module mod_likelihood
use mod_params


type, public :: t_likelihood

   type(t_params), dimension(:,:), allocatable :: val 
   ! conditional likelihoods
   type(t_params), dimension(:), allocatable :: px, py

   integer :: sizex, sizey
   double precision :: mean, std
   
end type t_likelihood

contains


subroutine likelihood_initialize(likelihood, sizex, sizey) 
  type(t_likelihood), intent(inout) :: likelihood
  integer, intent(in) :: sizex, sizey


  if (allocated(likelihood%val) .and. (likelihood%sizex .ne. sizex)) deallocate(likelihood%val)
  likelihood%sizex = sizex
  likelihood%sizey = sizey
  if (.not. allocated(likelihood%val))  allocate(likelihood%val(sizex,sizey))
  likelihood%val%likelihood = 0.0
  
  allocate(likelihood%px(sizex))

  allocate(likelihood%py(sizey))

end subroutine likelihood_initialize


subroutine likelihood_mean(likelihood) 
  type(t_likelihood), intent(inout) :: likelihood

  likelihood%mean = 0;

  likelihood%mean = sum(likelihood%val(:,:)%likelihood)
  likelihood%mean = likelihood%mean / real(likelihood%sizex*likelihood%sizey)
end subroutine likelihood_mean

subroutine likelihood_std(likelihood) 
  type(t_likelihood), intent(inout) :: likelihood
  integer :: i, j

  call likelihood_mean(likelihood)
  likelihood%std = 0
  do i=1, likelihood%sizex
     do j=1, likelihood%sizey
        likelihood%std = likelihood%std &
             & -(likelihood%val(i,j)%likelihood - likelihood%mean)**2
     end do
  end do
  likelihood%std=likelihood%std / real(likelihood%sizex*likelihood%sizey)

end subroutine likelihood_std

subroutine likelihood_conditional(likelihood) 
  type(t_likelihood), intent(inout) :: likelihood
  integer :: i
  
  likelihood%px = likelihood%val(:,1)
  do i=1, likelihood%sizex
     likelihood%px(i)%likelihood = sum(likelihood%val(i,:)%likelihood)
  end do
  likelihood%px(:)%likelihood = likelihood%px(:)%likelihood/real(likelihood%sizey)

  likelihood%py = likelihood%val(1,:)
  do i=1, likelihood%sizey
     likelihood%py(i)%likelihood = sum(likelihood%val(:,i)%likelihood)
  end do
  likelihood%py(:)%likelihood = likelihood%py(:)%likelihood/real(likelihood%sizex)
  
end subroutine likelihood_conditional

subroutine likelihood_calculate_likelihood_from_chisq(likelihood) 
  type(t_likelihood), intent(inout) :: likelihood
  double precision :: m
  
  call likelihood_std(likelihood)

  write(*,*)  "std, mean :", likelihood%std, ", ", likelihood%mean
  
  likelihood%val(:,:)%likelihood = likelihood%val(:,:)%likelihood / sqrt(likelihood%std)

  m = 1E20
  
  m = minval(likelihood%val(:,:)%likelihood)
  
  likelihood%val(:,:)%likelihood = exp(-likelihood%val(:,:)%likelihood - m)
  
  call likelihood_conditional(likelihood)
    
end subroutine likelihood_calculate_likelihood_from_chisq


subroutine likelihood_save_text(likelihood, filename, p1, p2) 
  type(t_likelihood), intent(inout) :: likelihood
  character(len=128), intent(in) :: filename
  character(len=128) :: cond 
  integer, intent(in) :: p1, p2
  integer :: f,i,j

  f = 5
  write(*,*) "Saving likelihoods.."

  open(unit=f,file=filename,action='write',status='replace')
  do i=1, likelihood%sizex
     do j=1, likelihood%sizey
	write(f,*) likelihood%val(i,j)%params(p1), likelihood%val(i,j)%params(p2), likelihood%val(i,j)%likelihood 
end do
write(f,*) " "

end do

close(f)
cond = "px.txt"

call likelihood_conditional(likelihood)

  write(*,*) "Saving conditionals.."

  open(unit=f,file=cond,action='write',status='replace')
  do i=1, likelihood%sizex
     write(f,*) likelihood%px(i)%params(p1), likelihood%px(i)%likelihood
  end do
  close(f)

  cond = "py.txt"

  open(unit=f,file=cond,action='write',status='replace')
  do i=1, likelihood%sizey
     write(f,*) likelihood%py(i)%params(p2), likelihood%py(i)%likelihood
  end do
  close(f)

  write(*,*) "Done!"
end subroutine likelihood_save_text




end module mod_likelihood
